(* orthogonal L^2 projection of vector c onto the L^1 unweighted ball 
   with radius tau : O(n*log(n)) version *)
let projector tau c =
  assert( tau >= 0.0 );
  (* test if c is already in the ball *)
  if ( Lacaml_ext.D.Vec.nrm1 c <= tau ) then
    Lacaml.D.copy c
  else
    let n = Lacaml.D.Vec.dim c in
    (* test if the weighted ball is reduced to 0 *)
    if ( tau = 0.0 ) then
      Lacaml.D.Vec.make0 n
    else
      let a = Lacaml_L1.D.Vec.abs c in
      Lacaml.D.Vec.sort ~decr:true a;
      let s = ref 0.0 in
      let phi = ref 0.0 in
      let idx = ref 0 in
      while ( !phi <= tau && !idx < n ) do
        idx := !idx + 1;
        s := !s +. a.{!idx};
        phi := !s -. (float_of_int !idx) *. a.{!idx};
      done;
      assert( !idx > 1 ); (* car on a éliminé le cas tau = 0.0 *)
      let lambda =
        if ( !phi > tau ) then
          let theta  = ( !phi -. tau ) /. (float_of_int (!idx - 1) ) in
          a.{!idx} +. theta
        else
          begin
            assert( !idx = n );
            ( !s -. tau ) /. (float_of_int n )
          end
      in 
        Lacaml_L1.D.Vec.soft_threshold lambda c;;
      
